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ABSTRACT 



We present a new code for the calculation of the ID structure and synthetic spectra of accretion disks. The code is an extension of 
the general purpose stellar atmosphere code PHOENIX and is therefore capable of including extensive lists of atomic and molecular 
lines as well as dust in the calculations. We assume that the average viscosity can be represented by a critical Reynolds number in a 
geometrically thin disk and solve the structure and radiative transfer equations for a number of disk rings in the vertical direction. The 
combination of these rings provides the total disk structure and spectrum. Since the warm inner regions of protoplanetary disks show 
a rich molecular spectrum, they are well suited for a spectral analysis with our models. In this paper we test our code by comparing 
our models with high-resolution VLT CRIRES spectra of the T Tauri star GQ Lup. 
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^ ' 1. Introduction 

Gas and dust disks are common objects; they can be observed 
around a variety of objects such as very young stars (e. g. T Tauri 
and Herbig Ae/Be stars), evolved binaries (cataclysmic vari- 
ables), and even black holes. With the discovery of the first 
extrasolar planets over 10 years ago, the interest in protoplan- 
etary disks has increased. Disk properties such as density, tem- 
perature, and chemical composition effect the process of planet 
formation and therefore also the characteristics of the planets. 
For very young stars (ages of ~ 1-10 Myr), the accretion 
of disk matter onto the star is an internal source of energy for 
the disk. The mass of such a disk is typically a few percent 
of the stellar mass and mass accretion rates are of the order of 
10"^ - 10"'*^ Mo/yr. Furthermore, irradiation by the central star 
heats the outer layers of the disk which, in turn, radiate the im- 
pinging energy away. Therefore, the direct detection of proto- 
planetary disks provides the possibility to observe the earliest 
evolutionary stages of planet formation. 

Tremendous efforts have been made to model the struc- 
ture and more importantly for the com parison with obse r vation s 
radiative transfer in accretion disks. Kriz & Hubenvl (Il986h . 
IShaviv & Wehrsd (Il986|), |Huben^ (tl990D and others went be- 
yond the ve rtically averaged mode ls of IShakura & SvunvaevI 
( il973h and iLvnden-Bell & Pringld (11974 or models using 
the diffusion approximation (e. g. iMever & Mever-Hofmeisterl 
Il982t ICannizzoetai]|1982h to obtain full numerical solutions 
for the structure of and the radiative transfer in accretion disks. 
Since then, these models have reached a high degree of so- 
phistication (for an ove rview of protoplanetary disk models see 
iDullemond et alj|2007h . 

Even though instruments like the VLTI constitute a major 
improvement for the observation of spatially extended objects, 
most of our information about the physics of protoplanetary 
disks comes from spatially unresolved spectra. NIRSPEC ob- 



servations of proto planetary disks hav e shown the capabilities 



of IR spectroscopy dNaiita et 



trograph CRIRES dKaeufl et 



tay 
[aD 



2003). 
2004) 



With the ESO IR spec- 
such observations have 
become available to a larger community of astronomers. Dust 
continuum radiative transfer (RT) calculations (e. g. Wolf 200l!) 
can reproduce the overall structure of the dust disk out to a few 
hundred AU. The more abundant gas in the disk, however, can- 
not be predicted by these models. The warm inner disks (tem- 
peratures between 100 K to a few 1000 K) provide a special 
laboratory to study the gas structure because temperatures and 
densities are adequate to produce molecular spectral lines visi- 
ble in the near- to mid-infrared. High-resolution observations in 
combination with model spectra enable us to obtain kinematic 
information about the gas since line profiles are governed by the 
velocity field in the disk. Another interesting observation is the 
agreement of mean inner gas disk rad ii and orbital radii of short- 
period extrasolar planets (ICarril20(j7l) . To further investigate this 
and other phenomena, detailed gas and dust models of the warm 
inner disk regions are necessary. 

Our paper is structured as follows: In Section 2 we will ex- 
plain the construction of our model structure and synthetic spec- 
tra. We will concentrate on the basic concepts and highlight the 
implementation of the solution. In Section 3 we will present syn- 
thetic spectra for the disk of the T Tauri star GQ Lup and com- 
pare these with CRIRES infrared observations to show the po- 
tential of our ID disk model code. Section 4 will give a summary 
and an outlook of our work. 



2. Models 

We have developed a circumstellar disk radiative transfer code 
as an ext ension of the well-te sted model atmosphere package 
PHOENIX dHauschildt & Baronll 19991) . PHOENIX can handle very 
large atomic and molecular line lists and blanketing due to sev- 
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eral hundred million individual lines is treated in the direct opac- 
ity sampling method. Dust is included in the models presented 
here by treating condensate formation under the assumption of 
chemical equilibrium and phase-equilibrium for several hundred 
species. I. e. all dust monomers exceeding the local saturation 
pressure as defined by thermal equilibrium are allo wed to con- 
dense (this implies a supersaturation ratio 5 = 1; AUard et"an 
1200 IL cf. also .Helling e t al. 2008. for a comparison of different 
condensation treatments). 

The grain opacity is calculated from the most important re- 
fractory condensates using optical data for a total of 50 differ- 
ent species. Absorption and scattering are calculated in the Mie 
formalism, assuming a given particle size distribution for a mix- 
ture of pure spherical grains, following the general setup of the 
Dusty set of PHOENIX atmosphere models jAUard etalJl200lL 
a plot with relative abundances of the most important species 
in our disk models is shown in Fig. |9]). This equilibrium as- 
sumption is a good approximation in the inner optically thick 
layers of the well-mixed disk atmosphere. In the low-density 
outer layers, non-thermal effects are more important; including 
these effects in the DISK version of PHOENIX is planned for 
the future. In the later phases of disk evolution grain growth 
will become important and may cause departures from a ho- 
mogeneous size distribution. Partial pressures for the different 
atoms, molecules, and dust species are calculated with the new 
"Astrophysical Chemical Equilibrium Solver" (ACES) equation 
of state (Barman, in preparation) which allows us to reach tem- 
perature regimes as low as a few tens of degrees Kelvin. 

In our models, the disk region considered is divided into 
rings (see Fig. [Hi and a plane-parallel disk atmosphere is com- 
puted between the midplane and the top of the disk for a given 
number of quadrature points (Gaussian angles) yu - cos 6 for 
each ring independently. Here 6 denotes the angle between the 
characteristic and the normal to the disk plane. 

We adopt the standard accretion model for geometrically 
thin disks, i. e. the disk height H is much smaller than the disk 
ring radius R. This assumption decouples the treatment of ver- 
tical and radial disk structure because the vertical structure is 
in quasi-static equilibrium compared to time scales for the ra- 
dial motion of gas. Matter is assumed to rotate with Keplerian 
velocities and viscous shear decelerates inner and accelerates 
outer parts leading to accretion of matter and outward trans- 
portation of angular momentum. Molecular viscosity is too small 
to provide the observed mass accretion rates. Thermal convec- 
tion in accretion disks was investigated by differ ent aut h ors us - 
ing various methods which are summarised in iKlahd (I2007h . 
The results show that thermal convection is unlikely the dom- 
inant source of turbulence and even leads to inward transport of 
angular momentum. Fur thermore, a heating source is required 
to drive the convection. iBell & LinI (1994) argue that convec- 
tive instabilities can only occur at temperatures T < 200 K or 
2000 K < r < 20000 K, i. e. temperature regimes which our 
models do not reach (see Sect.[3]l- The magnetorotational insta- 
bility (MRI) introduced by poloidal magnetic fields in weakly 
ionised disk matter (Balbus & Hawlev 1991) is the currently 
favoured origin of viscosity but its effect on the thermal struc- 
ture of the disk cannot be easily described or parametrised. Even 
though temperatures T > 1000 K are necessary to thermally 
ionise disk material, cosmic ray ionisation is possible at surface 
densities 2 < 100 g/cm^ (Gammie 1996) which is true for all 
our models presented below. The mean viscous dissipation is 
often modelled as an "alpha-viscosity" resulting in a vertically- 




Fig. 1. Disk ring structure as adopted for our calculations. The 
radius of the rings increases exponentially. The left panel shows 
a face-on view of a disk, the right one a vertical cut viewed edge- 
on (height is not to scale). The dotted lines are the radii R for 
which the models are calculated while the solid lines show the 
borders of a disk ring. The disk structure is assumed to be con- 
stant over the ring width. 



averaged viscosity 



(1) 



(I Shakura& SvunvaevllI973]) where < a < 1 is the angular 
momentum transfer efficiency, Cs the sound speed, and H the 
pressure scale height. Alternatively, a turbulent viscosity can be 
modelled assuming a critical Reynolds number Re 



yJGM^R 
R^ 



(2) 



(iLynden- Bell&Pringl^ll974 . This second model has the ad- 
vantage, that the calculation of the mean viscous dissipation is 
decoupled from the thermal structure of the disk and is adopted 
here. Both of these formalisms allow one to account for the effect 
of viscosity on the disk structure without the need to describe its 
origin in detail. 

2.1. Start models 

In order to start a new DISK disk calculation either an existing 
model can be used as input or a grey LTE start model is con- 
structe d. In the latter case we follow the approach of iHubenvl 
(119901) . 

The model is initialised by setting up a mass depth scale 
OT,, / = 1, . . . ,NL, where NL is the number of layers, which 
will be kept fixed also for later calculations with the same input 
parameters. The m scale is equally spaced on a logarithmic grid 
between the column mass at the midplane of the disk 



M 
6nv 



(3) 



and the outer value nii, which is an input parameter. This spacing 
is done for NL - Ni„ points, while for the remaining Nia layers 
each point j, j - 1 , . . . , Nin, is positioned halfway between Mq 
and nij-i, with mj^o = ^mnl-N;,,-!- This is done to provide nu- 
merical stability in the iterative process because of the steep m 
gradient near the midplane. A typical value for Nin is 6 or 12 for 
NL = 64. 

In a next step, we calculate the depth-dependent viscosity 



(m 



(4) 
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assuming a value of v determined by an assumed constant crit- 
ical^ Reynolds number (Eq. |2|, where ^ > is a free parameter 
jKriz & Hubenv 1986). The smaller the value of the lower the 
temperature. However, if irradiation is noticeably strong, the ef- 
fect o f 4" on the structure has neglig ible influence on the spec- 
trum. iLvnden-Bell & Pringld (1 19741) argue that the Reynolds 
number has to be chosen to equal the critical one for the on- 
set of turbulence. We usually take Re - 10"* - 10^ which leads to 
a = 0.1-0.01. 

We further assume an isothermal density structure - the 
sound speed associated with the gas pressure Cg and the flux 
mean opacity kh are independent of height - which lets one de- 
rive a simple analytic expression for the density at each depth 
point q(z). The relation 



m{z) 



r 



g{z) dz 



(5) 



is then inverted to convert our mass depth variable m into a 
height above the midplane of the disk z- 

Finally we employ the formal LTE solution derived 
in iHubenvl (.1 990) and the Rosseland opacity tables of 
[Ferguson et al.l (|20051) to get a temperature structure by iterating 



y^' 



eff 



T 1- 



2t,( 



V3 3eTtot V 



(6) 



with Kn = Kj^{T,g) until the T-r-structure is consistent for each 
layer. In Eq.|6]we set e = kr/kb equal to 1 and Ttot is th e optical 
thickn ess at the midplane of the disk (see Sec. Ill b of lHubenvl 
fT990l) . 



2.2. Iterative procedure 

After a restart model is read in or a start structure has been com- 
puted, the hydrostatic equation, the radiative transfer (RT) equa- 
tion, and the energy balance equation have to be solved. This is 
done iteratively until convergence in flux is obtained. Our termi- 
nation criterion is 



max 

= 1,...,NL 



Fmimd 



< IQ- 



(7) 



where F^ is the expected "mechanical" flux released by the 
disk and F^ the current flux value. This convergence criterion 
in Eq. [T] is not always reached and for a sufficiently small 
niax/=i .jvjL |Ar,| <K 1 K the calculation is then stopped after 
a given number of iterations. Typical errors in F are of the order 
of IQ-^. 



2.2.1. Hydrostatic equation 

The hydrostatic equation is one of the basic equations that gov- 
ern the structure of a stable disk. Since the mass of the disk is 
much smaller than that of the central object, we can neglect self- 
gravitation of the disk. Assuming that the radial component of 
the central star's gravitation and the centrifugal forces of the ro- 
tating disk just cancel each other out, we obtain the vertical hy- 
drostatic equation for a thin disk 

dP GM. 

where P is the sum of gas pressure f g and radiation pressure 
Pf. It is now convenient to use the relation dz/dm - -l/g and 
replace Eq. [8]by 

d^P _ cl GM^ 
dnfi 



P R' 



(9) 



This way, we eliminate the height z (which is not known a-priori) 
and introduce the sound speed = P/g, which is taken from 
the previous iteration and only depends on the temperature T. 
Therefore, the error in P and g cancel out and we obtain a more 
stable version of Eq.[8] This second-order differential equation is 
then decomposed into two coupled first-order differential equa- 
tions. With the inner boundary condition P'(Mo) - (symmetry 
at the midplane) a nd the outer boundary condition P(m\) derived 
in lHubenvl (Il990b . we solve the two point boundary value prob- 
lem with a simple shooting method. 

2.2.2. Radiative transfer 

We solve the plane-parallel radiative transfer equation 



dly 

H— ^ ly-Sy 

dTy 



(10) 



for a given number of Gaussian angles yu,-, / = 1 , . . . , NG, with 
corresponding quadrature weights a, from the midplane (z - 0, 
m - Mo) to the top of the disk (z - Zmax. = '^i)- The upper 
boundary condition for Eq. [TOlis 

Iy(v, Zmax) - /rV, -yU, 2max) (1 1) 

and the lower boundary condition due to symmetry conditions is 

/,(v,-;u,0) = «v,yu,0) . (12) 

In Eq. [TTl the expression /™'(v, -yu, Zmax) denotes the impinging 
radiation from the central star at the top of the disk atmosphere. 
The radiation field Jy is c omputed using the Accelera ted Lambda 
Iteration (ALI; see, e. g.. iHauschildt & Baronlll999l) . Hence, we 
have to incorporate the boundary conditions given in Eqs. [TT] 
and[T2]into the formal solution (FS) of the transfer equation ( fTOl i 
and into the Approximate Lambda Operator (ALO) A*, which 
is introduced to improve the convergence rate of the Lambda 
Iteration by reducing the eigenvalues of the amplification matrix. 
The FS can be written as 



>/new — AiSolds new — (1 ^)>^new + ^B , 



(13) 



where S is the source function and e is the thermal coupling 
parameter. The A-operator is split according to 



A = A* -H (A - A*) 

and we can rewrite Eq.[T3]as 

-/new - A*5 new + (A - h*)S old 



(14) 



(15) 



The non-local A*-o perator is calculated according to 
IHauschildt etaTI (11994 . 



2.2.3. Energy balance 

The standard accretion disk model demands that all the energy 
that is produced by viscous dissipation, i. e. mechanical energy, 
is released in form of radiative and convective energy, viz. 

£m(m) = E,{m) + Edm) . (16) 

For a viscous disk, the left hand term can be written as 



9 GM* 
Emim) = - -^-v{m)Q(m) 



(17) 



dKriz & Hubenvl[T986l) and the radiative energy is expressed by 

E^(m)^An \ (i]y(m) -Xy(m)Jy{m))dv (18) 
Jo 
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Fig. 2. Temperature correction scheme for a not irradiated disk atmosphere with R - 0.069 AU, M = 5 ■ 10"^ Mo/yr, Re - 5 ■ 10'^, 
R* - 1.57 Rq, and = 0.6 Mq. The top plot of the left panel shows the relative temperature change for each iteration in the top 
most layer The bottom plot is the temperature in that layer The right plot shows relative temperature corrections as a function of 
layer number for a chosen set of iteration steps (see legend). The dotted line shows the correction in the first and the dashed line in 
the last iteration. The correction is limited to AT/4 in the first four iterations to avoid overcorrections. In this example convergence 
in the Unsold-Lucy temperature correction is reached after ~ 35 iterations. 



We shall neglect the convective energy term Ec from here on 
because radiation is dominating the temperature structure and 
convection seems to have little effect dD'Alessio et a ljl998l) . We 
then derive an Unsold-Lucy class temperature correction scheme 
jLucv 1964) by integrating the second frequency-integrated mo- 
ment of the radiative transfer equation 



dK 

— ^khH 
am 

where 



/-*CO 

{XvlQ)HydvlH 

Jo 



The introduction of the Eddington factors 



Mm) = 



K(m) 
J(m) 



and fn = 



J(0) 



(19) 



(20) 



(21) 



and insertion of the first frequency-integrated moment of the ra- 
diative transfer equation gives the temperature correction law 



with 



4o-r3 



Kj ^ AH(O)fKiO) ^ 



Kb 



Ik/h 



JK Jo 



KHAH{m')dm' 



1 dAH(m) 



Kb dm 



(22) 



Kj 



/-*oo 

Jo 



/g) Jv dv/ J and kb 



/-*CO 

Jo 



/g)Bydv/B, (23) 



the absorption mean and Planck mean opacity. An example iter- 
ation history is shown in Fig.|2l 

2.2.4. Irradiation 

Irradiation by the central star plays an important role in the de- 
termination of the temperature and height profile of a protoplan- 
etary accretion disk. Therefore, we have taken special care to 
treat this effect in detail. The impinging intensity onto the sur- 
face of the disk (see Eq. [TTI ) is determined by first calculating 



the slope of the disk surface cos(^) — Az/AR at ring radius R 
by computing the height of the disk according to our start model 
(see SectionIO at R' ^ R + AR with AR ^ R ■ 10 -. For each 
Gaussian angle /i, and corresponding integration weight a,- the 
projected surface fraction on the star A, is determined. This area 
is then subdivided and a mean intensity reaching the disk surface 
is evaluated for a star with a blackbody spectrum according to 



/v(jUi, V) = — 

2n 



( NS 

\]=\ ^ ■' 



RA^ lyi^ipV)^ 



/v(l,v) 



with A,- = ^ 



(24) 



where NS is the number of surface segments for Gauss angle ju,, 
rj is the distance from the area element Aj to the disk surface, 
and the angle under which the radiation leaves the star. The 
factor 1 /27r is needed since the the irradiation onto the disk is 
unidirectional and not isotropic. Limb darkening is considered 
by the factor Iy{pLj,v)IIy(l,v). Alternatively, a PHOENIX spec- 
trum can be used as irradiation source. The irradiation geometry 
is shown in Fig. [3] 

2.2.5. Spectrum 

Having determined a set of self-consistent disk ring structures, 
we calculate a last iteration with a larger number of Gaussian 
angles (usually NG = 16) and a fine wavelength spacing (up 
to AA - 0.2 A in the wavelength range of interest) to obtain 
a high-resolution model spectrum. Since the Keplerian rotation 
velocity and the inclination angle / of the disk are not considered 
in the single ring spectra yet, we compute a combined spectrum 
according to 



J«'" Jo 



/v(v, = COS(0 



, r, i) r d(p dr 



(25) 



NR r 

.cos(o2 [(^rf -K")' 



NA 



Y^l[v'(i,(l)k),<l>k,Rj,i) 



k=l 



(26) 
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0,06 0,05 0.04 0,03 0.02 0,01 4,7265 4.7266 4.7267 4,7258 

z/R A [/jm] log(m/[g cm ]) 



Fig. 4. Optical depth structure for a disk ring model with R = 0.065 AU, = 0.8 Mq, = 2.55 TJq, M = 3 ■ lO"*^ Mq yr"\ and 
/?e = 5 ■ 10"*. The left plot shows the run of Text with height (z/R = is the midplane of the disk). The asterisks denote the optical 
depth at the center of the line, the crosses are for a continuum point (see right panel). The dotted line at r = 1 marks the region 
where the line and the continuum become optically thick, i. e. where the radiation that we see comes from. The plot on the right is 
a temperature-log m-diagram showing at which temperatures and column masses the optical depth for the line (dashed line) and the 
continuum (dotted line) become unity. 



In Eq.|26]we assumed that the disk is axis-symmetric and that the 
intensity is constant for all radii between inner and outer radii for 
each ring. In addition to the integration over all disk rings, the 
influence of the disk's rotation on the line profile is taken into 
account by applying the Doppler shift 



v'(i, 0) = V 



1 - u(i, (f>)/c 



^ 1 4- u(i, 0)/c 



(27) 



to the line. Here the velocity u(i, (p) = MKepier sin(/) cos(0) is de- 
termined for a given inclination / and for a set of disk ring seg- 
ments with azimuthal angle 0. We use NA =100 disk ring seg- 
ments, i. e. 100 steps in 0, to determine the rotationally broad- 
ened line profile. This method is a simple way to determine line 
profiles for rotating accretion disks if t he lines originate in the 
very upper layers of the disk. However. iHorne &Marshl(ll986h 
noted that an anisotropy in the local emission pattern changes 
the global line shape: line photons are trapped in optically thick 
emission layers and can more easily escape in directions where 
there are larger Doppler gradients. The true consequence of this 
on the line shape will be discussed in a future paper (Hiigelmeyer 




Fig. 3. Irradiation geometry as adopted for our calculations. We 
consider a star with radius at a distance R from a disk ring 
with a height Zmux and a slope of the disk surface cos(i^). For 
each Gaussian angle - a + if and its corresponding integration 
weight fl, the projected surface fraction on the star is determined 
and the irradiation intensity is calculated. The angle under which 
radiation leaves the surface of the star is fij = a + 6. 



et al. in preparation) where we will present full 3D radiative 
transfer calculations in rotating accretion disks. 




4,72 4,74 4,76 

A + offset [/im] 



ring weigtifs 



Fig. 5. The left panel shows normaUsed disk ring spectra. 
Intensities and wavelength are offset for clarity. The right panel 
depicts bars corresponding in height to the contribution of each 
disk ring spectrum to the total spectrum. The bar for the stellar 
contribution has to be multiplied by a factor of 25. The spectra 
and the weights are grey-scale coded corresponding to their ring 
radius R. 



3. Synthetic spectra for GQ Lup 

We retrieved spectra of T Tauri stars taken with the high- 
resolution infrared spectrograph CRIRES at the VLT fro m the 
ESO Science Archive Facility (see'Pontop pidan et al.l2008L for a 
description of the observations). The observations were reduced 
using a combination of the CRIRES pipeline and our own IDL 
routines. The telluric absorption lines in the spectrum were re- 
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0.031 AU 0.137 AU 

0.045 AU 0-200 AU 

0.065 AU 0.290 AU- 

0.094 AU 0.422 AU 




0.0100 0.1000 1.0000 
column moss density [g/cm^] 



Fig. 6. Temperature structures for the eight disk rings calcu- 
lated for the spectral analyses. The temperature profiles marked 
with symbols, e. g. for disk rings with R - 0.031 AU and 
R = 0.045 AU, are not included in the best fit model. 



moved by using a telluric model spectrum (Seifahrt, in prepa- 
ration). From the spectra taken between April 22 and April 26 
2007, we selected the observation of the classical T Tauri star 
GQ Lup to demonstrate the applicability of our models to obser- 
vations. 



3.1. GQLup 

GQ Lup is a classical T Tauri star (CTTS) which is 
mostly kno wn for its recently disco vered sub-stellar companion 
GQ Lup B (iNeuhauser et al.l [20051) . The activity due to ongo- 
ing accretion of the CTTS makes a well constrained determi- 
nation of its physical parameters difficult. This becomes obvi- 
ous regarding the diff'erences in visual brightness of more than 
2 mag (Vmax - 11-33 mag and Vmin - 13.36 mag; Janson et"al] 
120061) . iBroeg et al.1 ( 120071) and ISeperuelo Duarte et all ( l2008l) 
(BR07 and SD08 from here on) both obtained photometric and 
spectroscopic data to determine rotational periods and vsin/ 
for GQ Lup A and also to derive other stellar parameters 
which we need as input for our models. Both authors assume 
an effective temperature fo r a K7 V star of Tetf = 4060 K 
dKenvon & Hartmannlll995h . 

SD08 obtained a radius of 7?* = 1.8 + 0.3 Rq for the K7 V 
star assuming a mean distance of 150 + 20 pc. From evolutionary 
tracks they derive a mass of M* - 0.8 ± 0.2 Mq. Together with 
their rotational period of 10.7 + 1.6 d and a spectroscopically 
determined vsin/ = 6.5 ± 2.0 km s ' they find an inclination 
angle of/ = 51° + 13°. 

BR07 measure a shorter photometric period of 8.45 + 0.2 d 
and a larger radius of - 2.55 + 0.41 Rq from T^ff and a 
luminosity L from evolutionary tracks. With these values and 



vsin; = 
/ = 27° 



6.8 

:5°. 



0.4 km s they predict an inclination angle of 



3.2. Models 

We calculated small model grids for the input parameters from 
the publications of BR07 and SD08 (see Sec. 13. Il l, i. e. using 
Teff = 4060 K, = 0.8 Mq, and radii of R^ = 2.55 Rq 
and R-t. = 1.80 Rq. The different radii lead to luminosities 
of L = 1.58 Lq and L - 0.79 Lq, respectively. We varied 



the mass accretion rate (M = 8 ■ lO"'" - 2 ■ 10"^ Mq yr"^) 
and the Reynolds number (Re - 5 -10"^ - 1-10^). For 
each parameter set we computed eight disk rings with R = 

0. 031, 0.045, 0.065, 0.094, 0.137, 0.200, 0.290, and0.422AU. 
The ring radius R increases exponentially to achieve a small 
temperature gradient from ring to ring. Figure |6] shows typical 
temperature structures for a set of disk rings. The decreasing 
temperature in the outer layers stems from the fact that our 
temperature is set by the gas which has a relatively low opacity 
at the wavelengths where stellar irradiation is the strongest. This 
effect becomes small at radii larger than R - 0.065 AU. 

Since we consider a disk that has recently formed from a 
protostellar cloud and yet not seen significant grain growth, we 
assumed a standard ISM type power law grain size distribution 
with base size - 6.25 ■ 10'^ ^um and an exp onent of -3.5. The 
abundances follow iGrevesse & Noelsl (Il993h . 

In Fig. |5] we show the normalised spectra of each disk 
ring and how much each spectrum contributes to the total disk 
spectrum in the wavelength region considered. From these ring 
weights one can see that the outer disk radius /?out contributes 
only very little to the ring integrated spectrum and that the ex- 
tension to even larger disk radii would have almost no measur- 
able influence. Furthermore, it becomes obvious that the inner 
rings, even though their surface area is much smaller than those 
of rings at larger radii, have a much higher weight. This is sim- 
ply because the irradiated and viscously heated atmosphere is 
much warmer closer to the star and therefore has a higher flux 
level than that of rings further out. Furthermore, the continuum 
becomes optically thin for disk parts with R ^ 0.200 AU which 
further decreases the flux level. The influence of the central star 
GQ Lup A on the total spectrum is accounted for by a black- 
body spectrum of 4060 K. The ratio of stellar and disk flux 
is /^star/^disk ~ 5 depending on the input parameters for the 
disk. Figure |5] also shows that CO emission lines disappear at 
R > 0.290 AU which corresponds to temperatures in the outer 
line emitting layers ofT< 500 K. 

3.3. Spectral fit 

The comparison of our disk model spectra to the CRIRES spec- 
trum of GQ Lup in the wavelength range A - 4.59 - 4.81 fim 
indicates that the inclination of / = 51° + 13° derived by SD08 is 
too high: the emission lines are much broader in the model than 
in the observation. Even if we set to large values > 0.200 AU, 

1. e. excluding the inner rings which have broader lines due to 
larger Kepler velocities, the line widths cannot be acceptably fit- 
ted. For the input parameters of BR07, we find a reasonably good 
fit to the CRIRES data (see Fig.[8]l. The line widths are best fit- 
ted with an inclination angle of / = 22°, which is the lowest 
possible value given by the errors in BR07. The Une strengths 
are best reproduced if we set the inner radius to Rin = 0.052 AU 
and the outer radius to Rqm = 0.500 AU for a mass accretion 
rate of M = 3 ■ lO"** Mq yr"' and Re - 5 ■ lO'* coiTespond- 
ing to ff » 0.05. The line wings of the fundamental transition 
lines (v = 1-0) close to the continuum are somewhat broader in 
the observation than in the model. This argues for contributions 
of disk regions with high Kepler velocities, i. e. smaller radii R. 
However, if we include disk rings with smaller radii, the temper- 
ature sensitive v - 2 - I CO emission lines are too strong in the 
model. Furthermore, the low J R- and P-branch CO v = 1-0 
lines around 4.66 fim are overestimated by the model, or put in 
other words, the increment of the line strengths is too strong in 
the model. For our best model fit, we use a model which fits the 
strength of the higher order CO v = 1 - lines because these are 
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Fig. 8. CRIRES spectrum of GQ Lup (grey line) with our best fit disk model (black line). The model is calculated for a disk between 
Risi = 0.052 AU and Ro^t - 0.500 AU with a niass a ccretion rate of M = 3 ■ 10"^ yr ' and a Reynolds number of Re - 5 ■ 10'^. 
The stellar parameters are those of Broeg et alJ ('2007): = 2.55 Rq and T^g - 4060 K. For the stellar mass we adopted the value 
of Mi, - 0.8 Mq given by lSeperuelo Duarte et al . ( 2008). The hydrogen Pf/3 (7-5) line at 4.654 fim cannot be attributed to the disk 
because its Doppler width corresponds to a Kepler velocity of R = 0.025 AU - 2 Ri, which lies well within our determined inner 
disk radius. 



more frequent in the observation accepting a less well reproduc- 
tion of the low order lines. 

The part of the disk with radii R > 0.500 AU has no mea- 
surable influence on the spectrum at wavelength < 5 ;um. The 
strong CO fundamental transition lines (v = 1 - 0) in the spec- 
trum are well reproduced by the model. The weaker v = 2 - 1 
CO emission lines are also visible in the spectrum and fitted by 
our model. The '^CO isotope is weakly present in the model, 
which assumes the solar system carbon isotope ratio of 89:1, 
but we cannot find clear evidence for its existence in the ob- 
served spectrum. Woods & Willacv (2009) have modelled iso- 
topic fractionation for a protoplanetary disk quite similar to our 
GQ Lup model, but found only mild efifects on carbon monox- 
ide, with very modest increases in the '^CO/'^CO ratio only in 
the outermost layers of the disk. Considering the noise in the 
spectrum and the blending of '^CO and '^CO lines we there- 
fore cannot find evidence for a diff'erent isotope ratio in our ob- 
servations. The maximum temperature in the surface layers of 



the disk atmospheres is T = 950 K which is much smaller than 
the estimated excitation temperature of Tex = 1814 K found by 
Naiita et al. (2003) for DF Tau which has a similar spectrum but 
stronger v = 2 - 1 emission lines. 

In Fig.|2]we show the influence of the dust grain size dis- 
tribution on the spectrum. All of the spectra are based on the 
same structure model but the emergent flux is calculated for four 
different dust grain size setups, i. e. no dust, ISM grain size dis- 
tribution with gq - 6.25 ■ lO"-' /im, ao = 6.25 • 10"^ /zm, and 
flo = 6.25 ■ 10""^ yum. As we can see, neglecting dust leads to a 
lack of opacity and too strong lines. The ISM grain size distrbu- 
tion fits the data equally well as the smaller grains. The spectrum 
with ten times larger grains than ISM still reproduces the obser- 
vation reasonably well, but lines are slightly stronger than for 
ISM sizes. 

An interesting feature in the observed spectrum is the hydro- 
gen PfjS (7-5) emission line. This line has been attributed to an 
absorption line in the telluric standard star used for the calibra- 
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Fig. 7. Spectra calculated from the same structure model but for 
different dust grain sizes. Neglecting dust as opacity source re- 
sults in overestimation of the lines. ISM grain size distribution 
and a dust grain size setup with ten times smaller dust radii yield 
a very similar spectrum. The dust opacity is slightly reduced if 
we increase the grain size by a factor of ten. 
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Fig. 9. Relative abundances of the most important dust con- 
tributors to the opacity. The left plot is for the disk ring with 
R = 0.094 AU and the right one for R = 0.290 AU. Forsterite 
(Mg2Si04) is the most abundant dust species and strongly 
present in all layers. Graphite sets on around m = 0.01 g/cm^ 
where it blocks stellar irradiation efficiently and therefore sup- 
presses heating of deeper layers. 

tion of CTTS by 'Najita et alJ (l2003h . Since we used a telluric 
model to remove telluric absorption lines, we can relate the Pf/3 
feature unambiguously to an emission line intrinsic to GQ Lup. 
We measure a blue shift of 7 km/s for the line relative to the 
star and the width corresponds to a velocity of m a; 172 km/s. 
This value is much larger than the measured vsin/ = 6.8 km/s 
by BR07. Therefore, the line cannot be attributed to the star di- 
rectly. If we assume Keplerian rotation, the formal origin of the 
line is at 7? - 0.025 AU = 2 and could originate from the 
accretion flow onto the star or a stellar or disk wind. 

4. Summary & outlook 

We present a new ID structure and radiative transfer program 
for circumstellar disks that is an extension of the general purpose 



stellar atmosphere code PHOENIX. The disk is separated into con- 
centric rings and the structure is calculated from the midplane 
of the disk to the top of the disk atmosphere for each ring. A 
combination of all rings then yields the total disk spectrum. Our 
program assumes a geometrically thin accretion disk geometry 
and a critical Reynolds number to describe the energy release in 
the disk due to turbulent viscosity. We have modified the hydro- 
static equation, the inner and outer boundary condition of the RT 
equation, as well as the energy balance equation to account for 
the difference between star and disk atmosphere. Irradiation of 
the central star onto the top of the disk atmosphere is taken into 
account. 

Our disk model spectra can reproduce high-resolution IR 
emission line spectra, in this case of the CTTS GQ Lup. For 
this particular object, we calculated a set of disk models for 
different stellar inpu t param eter s given by the recent publica- 
tions of Broeg et alj (l2007l) and ISeperuelo Duarte et al (2008|) 
and varied the mass accretion rate M, the Reynolds number Re, 
and the inner and outer disk radii. We investigated the contri- 
bution of each disk ring on the total disk spectrum and showed 
where line and continuum radiation originates in the atmosphere. 
Furthermore we could show that the inclination of the warm in- 
ner disk of GQ Lup is ~ 22° which is much s maller than the 
value obtained bv lSeperuelo Duarte et alj (l2008l) and jus t withi n 
the uncertainties given for the inclination bv'Br oeg et al.l ( l2007h . 

Even though we could provide a reasonable model fit to the 
observed spectra of GQ Lup with our ID structure and RT code, 
there is room for improvement. One obvious drawback of our 
code is the assumption of equilibrium chemistry and local ther- 
modynamic equilibrium. This is not very likely to be valid for 
the cool inner and optically thin and irradiated outer disk atmo- 
spheric layers. Therefore, the consideration of departure from 
chemical equilibrium due to turbulent mixing or photodissocia- 
tion, and radiative excitation of molecules needs to be considered 
in future projects. Furthermore, convective energy transport will 
be considered in the future in order to investigate the influence 
of this effect on the disk structure and spectra. 

We are extending our effort to model the structure and spec- 
tra of circumstellar disks to full 3D radiative transfer calcula- 
tions. This way we can relax the assumption that there is no flux 
exchange between matter at different radii R and we will "natu- 
rally" consider the direct irradiation of the central star onto the 
disk and heat the very inner disk rim. The influence of the Kepler 
velocity field in the disk on the line profile will be investigated by 
means of two level atom line transfer (see iBaron & Hauschildn 

llool . 
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